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In this paper we present a hydrodynamic approach to de- 
scribe the motion of migrating bacteria as a special class of 
self-propelled systems. Analytical and numerical calculations 
has been performed to study the behavior of our model in 
the turbulent-like regime and to show that a phase transition 
occurs as a function of noise strength. Our results can ex- 
plain previous experimental observations as well as results of 
numerical simulations. 



I. INTRODUCTION 

Bacterial colony formation has been the subject of re- 
cent studies performed by both microbiologists Jj]-(j| and 
physicists |6H13|1 with the aim to acquire a deeper under- 
standing of the collective behavior of unicellular organ- 
isms. 

Since bacterial colonies can consist of up to 10 10 mi- 
croorganisms, one can distinguish several length and time 
scales when describing their development. At the largest, 
macroscopic scale (above cca. 10 -3 m) one is concerned 
about the morphology of the whole colony, which can 
be often described in terms of fractal geometry . 
At this length scale, the physical laws of diffusion gov- 
ern colony formation, and only certain details of the un- 
derlying microscopic dynamics (such as anisotropy Jl4| ) 
can be amplified up to this scale by various instabilities. 
Thus these large-scale features could have been success- 
fully understood using relatively simple models which ne- 
glect most of the microscopic details. In contrast, when 
looking at the microscopic length scale of individual bac- 
teria (cca. 10" 6 m), such global physical constraints play 
a much smaller role and it is the biology which provides 
us the knowledge of how the microorganisms move, mul- 
tiplicate and communicate. In this work we focus on the 
intermediate length and time scales which range from 
10~ 5 m to 10~ 3 m and from 1 s to 10 2 s. In this regime 
the proliferation of the bacteria can be neglected, and 
the most striking observable phenomena are connected 
to the motion of the organisms. 

Bacterial swimming has been well understood in liquid 
cultures [[D5|~[r7|, where the interaction between cells is 
usually negligible and the trajectory of each organism can 
be described as a biased random walk towards the higher 
concentration of the chemoattractants, e.g., certain nu- 
trients. In fact, the bacterial motion consists of straight 
sections separated by short periods of random rotation 



called tumbling, when the flagellas operate in a quite 
uncorrelated manner. If an increase of the local con- 
centration of the chemoattractant is detected, bacteria 
delay tumbling and continue to swim in the same direc- 
tion. This behavior yields the above mentioned bias to- 
wards the chemoattractant. Recent experimental obser- 
vations |3],|Il],^j revealed that under stress conditions even 
chemotactic communication occurs in colonies: chemoat- 
tractants or chemorepellents are emitted by the bacte- 
ria themself, thus control ling the development of various 
self-organized structures jL2[ . 

However, numerous bacterial strains like Pro- 
teus mirahilis, Serratia marcenses ^ Bacillus circu- 
lans, Archangium violanceum, Chondromyce s ap iculatus, 
Clostridium tetani [ [l8| or Bacillus subtilis |l0f are also 
able to migrate on surfaces to ensure their survival un- 
der hostile environmental conditions. In all of these ex- 
amples, irrespective of the cell membrane (Gram + or 
— ) and the type of locomotion (swarming or gliding), 
a significantly coordinated motion can be observed: the 
randomness, which is a characteristic feature of the swim- 
ming of an individual organism, is rather repressed. In 
particular the swarmer cells of the strain Proteus are very 
elongated (up to 50/im) and move parallel to each other 
Q . The strain Bacillus circulans obtained its name from 
the very complex flow trajectories exhibited, which in- 
clude rotating droplets or rings made up of thousands 
of bacteria [[b|,[l9| ■ To some extent similar behavior can 
be observed in each of the above listed examples, which 
naturally raises the long standing question as to whether 
this form of migration requires some external (or self- 
generated) chemotactic signal or whether simple local 
cell-cell interactions are sufficient to explain such cor- 
related behavior. 

One of the purposes of this paper is to show that the 
answer to such questions requires taking into account 
that the bacteria are far-from equilibrium, self-propelled 
systems. Recent studies have revealed that such systems 
show interesting and unexpected features: in various traf- 
fic models spontaneous jams appear and there is 
a strong dependence on quenched disorder p^ ]. Many 
biological system (such as swarming bacteria, schools of 
fishes, birds, ants, etc.) can be described as special self- 
propelled systems with local velocity- velocity interaction 
introduced by Vicsek et al. [^3| and also studied in a 
lattice gas model HJ. Both of these models revealed 
transition from disordered to ordered phase. Although 
most of the models above incorporate some level of dis- 
creteness, a continuum approach also has been applied to 
study the traffic flow problem |21| or to describe the geo- 
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FIG. 1. Schematic crossection of a bacterial colony. 



tactic motion of bacteria |25| . Toner and Tu have shown 
using renormalization group theory that in a closely re- 
lated continuum model ]26[ ] a long range order can de- 
velop even in two dimensions; this is not possible in the 
equivalent equilibrium models. 

In this paper we construct a continuum description for 
the particularly interesting collective motion in bacterial 
colonies. First we introduce our model, then we give a 



theoretical analysis of the problem. In Section |V| we 
present the numerical results and finally in Sectionjy] we 
draw our conclusions. 



II. THE MODEL 

When migrating on surfaces, bacteria usually form a 
very thin film consisting only of a single or a few lay- 
ers of cells, thus their motion can be regarded as quasi 
two-dimensional. Due to the small length scale, the very 
small Reynolds number yields an overdamped dynamics 
in which the thermal fluctuations (Brownian noise) and 
the various surface effects (e.g., surface tension, wetting) 
dominate the motion. The constrains of the cell-cell in- 
teraction must also be taken into account since the elon- 
gated migrating bacteria usually align parallel to each 
other, which results in an effective viscosity that tends 
to homogenize the velocities as we show later. Bacteria 
also tend to move close to the substrate which manifests 
in an effective gravity, flattening the colony. 

The particles of the "bacterial fluid" (i.e., the colony) 
are bacteria surrounded by a droplet of moisture (ex- 
tracellular fluid, see Fig. |l|) which is essential for them 
both to move and to extract food from the substrate. 
The volume (V*) and the mass (m*) of such a particle is 
considered to be constant, which gives a constant (3D) 
density g* = m* /V* . A typical bacterial colony grown 
on a surface can be considered as flat, so our fluid will 
be characterized by the 2D number density g(r, t) and 
velocity v(r, t) fields, where r S R 2 and t is the time. 
The number density is simply related to the local height 
of the colony h by h(r, t) — g(r, t)V* . 

In contrast to fluid mechanics, here the fluid element 
consists of only a few "atoms" . Thus random fluctuations 
do not cancel out and they have to be taken into account 
by an additional noise term in the equation of motion for 



the velocity field. The two basic equations governing the 
dynamics are the continuity equation 



d t g + V{gv) =0 



(1) 



and the equation of motion, which is a suitable form of 
the Navier-Stokes equation 

<9 t v + (v • V)v = — -Vp + iA7 2 v + — F(v) - -v 

g* g* t 

+V(*,t), (2) 

where p is the pressure, v is the kinematic viscosity, F(v) 
is the intrinsic driving force of biological origin, r is the 
time scale associated with substrate friction and 7](pc, t) 
is an uncorrelated fluctuating force with zero mean and 
variance a representing the random fluctuations. Let us 
now go through the terms of Eq. |2|. 

The pressure is composed of the effective hydrostatic 
pressure, the capillary pressure, and the externally ap- 
plied pressure 



p = g*gh + 7 V 2 h + p cxt . 



(3) 



If the radius of surface curvature is larger then the cap- 
illary length 
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then the capillary pressure can be neglected. Since this 
is normally the case (l ca .p ~ 3 mm for water) , we consider 
only the first term of Eq. || in the rest of the paper. The 
surface tension becomes relevant only at the boundary of 
the colony where the local curvature is not negligible. 

The viscous term in Eq. || represents not only the vis- 
cosity of the extracellular fluid, but also incorporates the 
local ordering of the cells. As we show in Appendix A for 
a simple hard rod model of the cell-cell interaction, the 
deterministic part of the dynamics of the orientational 
ordering is described in a similar manner to p3|: 



d^ 
dt 



/' 



(5) 



where the orientation of the zth rod is denoted by < 
9i < 7r and (-) e denotes spatial averaging over a ball of 
radius e. The angle 9 can be readily replaced by eg, 
a unit vector at angle 9 to the x axis. If the changes 
in the magnitude of the velocity are small then eg can 
be replaced by v and Eq. ^ yields an interaction term 
proportional to (v) e — v. 

Taking Taylor series expansions for the velocity and 
the density fields yields 



4l <e ddve+(£V)v e +i(£V) 2 ve- 
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If the density changes are small, we recover the viscous 
term of Eq. ^, so it, in fact, includes the self- alignment 
rule of previous models. 

Bacteria tend to maintain their motion continuously 
by propelling with their flagcllas. This rather complex 
behavior can be taken into account as a constant magni- 
tude force acting in the direction of their velocity 



*c v 



(6) 



where c is the speed determined by the balance of the 
propulsion and friction forces. That is, c would be the 
speed of a homogeneous fluid. 

Combining Eqs. |, §, § and § we obtain the following 
final form for the equations of the bacterial flow: 



d t h + V{hw) =0, (7) 



and 



d t -v + (v • V)v = -gVh + -Vpoxt + vV 2 v 

Q 

c v 1 

+-T-T --V + J7- (8) 

T |V| T 

These equations are similar to those studied in [^6| but 
here we derived them from plausible assumptions based 
on the underlying microscopic dynamics instead of phc- 
nomenological concepts. Nevertheless, the analysis of 
Toner and Tu also holds for our equations which per- 
mits us to use our model as a test of their prediction of 
the existence of a phase transition. 



III. ANALYTICAL RESULTS 

For certain simple geometries of the boundary condi- 
tion it is possible to obtain analytical solutions for the 
noiseless (u = 0) stationary state of our model if we sup- 
pose incompressibility (h = const.). Taking dt = in 
Eqs. and g[ the following dimensionless equations are 
obtained: 



VV = 



(9) 



|v(r)| = 1, 



(11) 



where the orientation of v is arbitrary, but independent of 
the spatial position. In this state there is a sustained non- 
zero net flux of the fluid so it is regarded to be ordered. In 
the next section we examine numerically how increasing 
the noise level yields a disordered state where there is no 
net flux present. 

Another simple, but practically more relevant bound- 
ary condition is realized when the system is confined to 
a circular area of radius R. In contrary to the previous 
example this is a finite geometry which in turn means 
that in the stationary state no net flux is possible. Thus 
the flow field must include vortices and we show that a 
single vortex is indeed a possible stationary configuration 
of the system. 

Let us assume that the velocity field, expressed in polar 
coordinates (r, 4>), is a function only of r: 



v(r) 



(12) 



where is the tangential unit vector. Taking the expres- 
sion of the nabla operator in polar coordinates, one can 
easily see that Eq. ^ is satisfied for any velocity profile 
v(r). Substituting Eq. [l^ into Eq. 10 yields two ordinary 
differential equations: 



= r z 



,d*v 
dr 2 



r- v(l + r z ) + r l 

dr 



and 



1 dp cx < 
g*X dr 



(13) 



(14) 



The first equation gives the velocity profile while the sec- 
ond determines the pressure to be applied for maintaining 
the constant height (h =const.) condition. 

The boundary conditions for the velocity profile are 
either v(0) = and v(R') — (closed boundary at 
r = R' = R/X) or v(0) = and ^(R') = (free bound- 
ary at r = R'). The homogeneous solution of the equa- 
tion of motion with the above boundary conditions is 
given by Ii(r), the modified Bessel function of order one. 
The particular solution is — %Li(r), where L\ is the mod- 
ified Struve function of order one [£7J . Thus the velocity 
profile of the single vortex stationary state in a noiseless 
system with circular boundary is 



and 



CT 



A 



(v'V')v' 



Vp 



V'V 



(10) 



where v' = v/c, A = y/ur and the V operator derivates 
with respect to r' = r/A. For the sake of simplicity we 
drop the prime in the rest of this section. 

First let us consider the simplest geometry when the 
system is defined on an infinite plane. In this case the 
stationary state is trivially 



v(r) = ah(r) - ^L^r). 



(15) 



The parameter a should be chosen to satisfy the bound- 
ary condition at r = R' . In Fig. || we show velocity 
profiles for different values of A having R = 1 fixed and 
imposing closed boundary condition. The maximal ve- 
locity decreases for decreasing R/X ratio (Fig. ||); thus 
the system behaves more similar to the usual (not self- 
propelled) systems where v(r) = 0. In this respect R/X 
is a measure of the self-propulsion. From Fig. | it is also 
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FIG. 2. Velocity profiles in a vortex for various values of 
Re' (bottom curve: Re' = 1, top curve: Re' = 16). 

clear that the minimal size of a vortex is A so if R/X ^> 1 
then many vortices are likely to be present in the system. 
This phenomenon is analogous to the appearance of tur- 
bulent motion in fluid dynamics and thus R/X can be 
regarded as a quantity analogous to the Reynolds num- 
ber: 

Re' = ^. (16) 
X 

Another dimensionless quantity characterizing the flow 
is ct/A which gives the relative strength of the inertial 
forces compared the viscous ones. This quantity can be 
regarded as an internal Reynolds number 

Re" = 2L. ( i 7) 
In our simulations we always had Re" -C Re'. 

IV. NUMERICAL RESULTS 

For further investigations we used numerical solutions 
of Eq. fj] and Eq. ||. To study the closed circular geome- 
try we implemented our model on a hexagonal region of 
a triangular lattice with closed boundary conditions. We 
used a simple explicit integration scheme to solve numer- 
ically the equations. 

We started our simulations from a uniform density 
and a random velocity distribution. Figure ^| shows the 
stationary state for the high viscosity (A ~ 3.16) and 
high compressibility (g — 750) case, corresponding to 
Re' ~ 4.4. The length and direction of the arrows show 
the velocity, while the thickness is proportional to the 
local density of the fluid. In Fig. g we present the ra- 
dial velocity distribution for the vortex shown in F ig. H 
and the velocity profile given by our calculations (Eq. [l5^ T 



Re' 

FIG. 3. Position (solid line) and value (dashed line) of 
the maximal velocity for the profiles on Fig. |as a function 
of Re'. At the maximal velocity, open boundary condition 
(dv/dr — 0) is satisfied. 

Rather good agreement is seen; the differences are due to 
the fact that our numerical system is not perfectly circu- 
lar. We have also performed velocity profile analysis for 
two vortices from [[19[ . The data and the fitted curves are 
displayed in Fig. |6|. Again, reasonable agreement is seen 
which means that Eqs. ^ and |§| give a correct continuous 
description for this particular example of self-propelled 
systems. 

In the case of real bacteria one hardly can observe a 
perfect vortex, so we tuned the parameters to see the be- 
havior of the model far from stationarity. Lowering the 
value of the compressibility g, we observed temporally 
periodic structures instead of a constant velocity profile 
vortex. Figure ^ shows such a configuration. The den- 
sity at the lower left part of the system is higher than 
at the top right part. As the system evolves in time 
the whole flow pattern rotates anti-clockwise, thus this 
state has non-zero net flux which oscillates in time. Sim- 
ilar behavior also has been reported for certain bacterial 
colonies Q . 

Another interesting case occurs when the compressibil- 
ity is high and the viscosity is low, which corresponds to 
high values of Re . In this limit we observed a long life- 
time multi- vortex state (Fig. |^). For the parameter val- 
ues we use, these vortices eventually disappeared leaving 
behind a single vortex. We conjecture that there exists 
a Re c above which the multi-vortex state does not die 
out, this would correspond to the turbulent state of the 
normal hydrodynamics. 

We also performed simulations in the open planar ge- 
ometry. In this case our numerical system was confined 
to a square with periodic boundary conditions. Our goal 
was to find if there exists the phase transition as a func- 
tion of the noise strength a. To do this we first defined 
an order parameter 

$ = <|v|), (18) 
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FIG. 6. Velocity profile data for two vortices taken from 
[ pi[ (circles and squares) and fitted profiles using Eq. [l^. 



FIG. 4. A numerically generated vortex for Re' ~ 4.4 and 
g — 750. The length of the arrows is proportional to the local 
velocity, while their thickness is proportional to the density. 
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FIG. 5. The measured (circles) and the theoretical (solid 
line) velocity profile for the vortex in Fig. W. 



FIG. 7. Non-stationary configuration at Re' = 4.4, g 



-■ 300. 
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FIG. 8. Multiple vortex state at Re' = 25.3, g = 128. 

where (•) denotes spatial average over the entire system. 
If $ = there is no net current and the system is in 
disordered state. If $ > some level of order is present. 
We have performed long time (> 10 7 Monte-Carlo steps) 
runs for two different system sizes (L=24 and 48) to check 
the presence of the transition. Figure || shows the results 
obtained which strongly suggests that there is a transi- 
tion around a ~ 6.5. This value is by far not universal as 
it depends on the actual values of parameters used in the 
simulations. However, the extraction of the critical ex- 
ponents would require much larger computational efforts 
(i.e., larger system sizes). 



V. CONCLUSION 

We have investigated a continuous model for coopera- 
tive motion of bacteria. For simple geometries and spe- 
cific parameters we were able to obtain analytical results 
for the velocity profile of vortices often observable on 
intermediate length scales in various bacterial colonies. 
Our results are in quantitative agreement with biologi- 
cal observations Q and with previous simulations |jl9| . 
We also showed that there exists a transition p(| from 
an ordered to a disordered phase as a function of noise 
strength. 

In its present form our model does not give a model for 
the pattern formation of bacterial colonies. However, it 
is rather straightforward to include surface tension and 
wetting effects in their full detail [£8| to handling the 
dynamics of the colony border. Including proliferation, 
the presented model can be also expanded towards larger 
time and length scales to study the fingering instabili- 
ties. It is also possible to introduce chemotactic response 
which would allow long range interactions to build up. 

Further work needs to be done to characterize the tur- 
bulent regime (which is mostly observed in swarming 
colonies) and to relate the effective parameters used in 
our model to the real world quantities such as food or 
agar concentration. 
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FIG. 9. The order parameter $ as a function of the noise 
strength a. 
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APPENDIX A: SIMPLE GEOMETRICAL MODEL 
FOR THE ORIENTATION AL ORDERING OF 
CELLS 



B 



In this appendix we show that the simple geometrical 
constrains originated from the rod-like shape of bacteria 
can motivate the local orientational ordering introduced 
in Sec. II. Microscopic observations of high density cell 
cultures reveal that the orientational order emerges since 
bacteria tend not to "intersect" each other: usually bac- 
teria form well defined layers in which they are parallel 
aligned. To investigate quantitatively this phenomenon 
let us consider a system of rods with unit length and neg- 
ligible width, covering the plane with a uniform average 
density g. Here we focus on the local orientational order 
only and neglect the translational movement of the cells. 
Thus let us assume that the only interaction in the sys- 
tem is the intersection of the rods, and the overdamped 
dynamics of the system is determined by two effects, a 
relaxation process minimizing the total number of inter- 
sections U and a random Brownian fluctuation as 



Ml A 9 11(9 9 



.) + mod 7r, 



(Al) 



where the direction of the i-th rod is denoted by 9i and 
£j is an uncorrelated noise: (£i(t)£j(t')) = CSijS(t — t'). 
As the rods are not directed, < 9 t < n. The quantity 
U can be written as a sum of "pair potentials" Vy : 



U(6 1 ,6 2 ,..) = lY, V v^ ei >*j> e ^ 



(A2) 



where Vy equals 1 if the rods i and j intersect and 
otherwise. 

As a mean-field approximation, we introduce the prob- 
ability density function P(0), which gives the probability 
density of the event that the direction of a randomly 
selected rod is in the interval [9,9 + d9}. We are inter- 
ested in the local (short range) ordering; thus the spatial 
homogeneity of P is assumed. Now 11(6), the expected 
number of intersections of a rod directed at angle 9, can 
be expressed as 



U(6) 



dV / d6'V(xi,6,x',6')P(6'), (A3) 



since gP(6) gives the probability density for the existence 
of a rod pointing in the direction 9 at any position. Due 
to the translational and rotational invariance of the sys- 
tem, U (9) is given in the form of 



U(0) = gf dd'P(6') I d 2 x'T/(0,0,x',6>' 
Jo Jv 



9). (A4) 



As J d 2 x.'V(0, 0, x', 9' — 9) gives the area of a parallelo- 
gram shown in Fig. ^ (g§ can be written as 



U(6) = g d6'P(6')\sm(0' 
Jo 



(A5) 



A 



FIG. 10. The rods A and B with a given orientation inter- 
sect each other if the center of rod B is inside of the displayed 
parallelogram attached to rod A. 



The time evolution of the probability density function 
i s g iven by the Fokker-Planck equation obtained from 
& 



dP 8 / B C 8 
__ | \ p TT a P 

8t 80 \ 86 2 89 



(A6) 



The steady state solution of Eq. A6 assuming detailed 
balance satisfies 



S 8U C8P n 
AF W + 2W = °< 



yielding 



P(6) = P exp 



2XU(0) 
C 



(A7) 



(A8) 



where the normalization for P(6) is provided by the pref- 
actor Po- 

Let us measure 9 compared to the spontaneously se- 
lected orientation (9). For small deviations U(9) can be 
expanded in a Taylor series as 



U(9)tt^6 2 + 0(6 4 ). 



(A9) 



Now, from Eqs. A8 and ^9 the coefficient A can be cal- 



culated for given g, A and D. Substituting Eq. AE into 
Eq. A8 gives 



P(6) w P e- Axe ^ c , 



with 



Po = 



(A10) 



(All) 



For 6, 9' < 1 the | sin(0' - 6)\ term of Eq. |A5| can be 
linearized both in 9 and 9': 



U(6) 




-7T/2 



d0'P(6') 



tt/2 



d6'P(6') 



d9'P( 



-7T/2 



7T/2 



d9'P(9' 



(A12) 
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After some algebraic manipulations the above expression 
can be simplified to 

r e 

U{9) = 2q d9'P(9')(9 - 9') » g P Q 9 2 (A13) 
Jo 

yielding the consistency condition gP = A/2. Thus, in 
this mean field approximation 



(9i - (0)i 



(A14) 



d^ _ _A_ 

dt ~~ttC 

that is we recover Eq. |5[ The angle 9 can be readily 
replaced hy eg, a unit vector at angle 9 to the x axis. If 
the changes in the magnitude of the velocity are small 
then eg can be replaced by v and we recover the self- 
aligning interaction Eq. [|. 
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